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Abstract 

In 1983 Felderhof, Ford and Cohen gave microscopic explanation of the famous Clausius-Mossotti 
formula for the dielectric constant of nonpolar dielectric. They based their considerations on 
the cluster expansion of the dielectric constant, which relates this macroscopic property with the 
microscopic characteristics of the system. 

In this article, we analyze the cluster expansion of Felderhof, Ford and Cohen by performing 
its resummation (renormalization). Our analysis leads to the ring expansion for the macroscopic 
characteristic of the system, which is an expression alternative to the cluster expansion. Using sim¬ 
ilarity of structures of the cluster expansion and the ring expansion, we generalize (renormalize) the 
Clausius-Mossotti approximation. We apply our renormalized Clausius-Mossotti approximation to 
the case of the short-time transport properties of suspensions, calculating the effective viscosity 
and the hydrodynamic function with the translational self-diffusion and the collective diffusion 
coefficient. We perform calculations for monodisperse hard-sphere suspensions in equilibrium with 
volume fraction up to 45%. To assess the renormalized Clausius-Mossotti approximation, it is 
compared with numerical simulations and the Beenakker-Mazur method. The results of our renor¬ 
malized Clausius-Mossotti approximation lead to comparable or much less error (with respect to the 
numerical simulations), than the Beenakker-Mazur method for the volume fractions below f ~ 30% 
(apart from a small range of wave vectors in hydrodynamic function). For volume fractions above 
f ~ 30%, the Beenakker-Mazur method gives in most cases lower error, than the renormalized 
Clausius-Mossotti approximation. 
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I. INTRODUCTION 


Einstein was the first, who applied statistical physics to calcnlate the viscosity of snspen- 
sion [1], Having in mind nanometer size sngar molecnles immersed in water, he considered 
a model of snfficiently big spherical particles immersed in viscons liqnid. Experiments show, 
that in this case the observed viscosity increases [2], In his work, Einstein related the 
observed (effective) viscosity of snspension with its structnre on the microscopic level. 
His resnlt, ?7eff/h = 1 + 5/20 - where rj denotes the viscosity of solvent and 0 denotes 
the volume fraction of the system - is valid only for dilute suspensions. This limitation 
is caused by the assumption, that the particles immersed in fluid do not influence their 
mutual motion. The problem of the influence of the particles on their mutual motion is 
essential to go beyond the diluted regime and was already addressed by Smoluchowski. 

His analysis for two sedimenting spheres leads to the follow¬ 
ing conclusions. The two spheres sediment faster than a single 
one. Moreover, the velocities of both spheres are deviated from 
the direction of the gravity held, as shown in the hgure. This 
example clearly demonstrates, that two sedimenting particles 
in gravity held behave diherently than a single particle, be¬ 
cause a single particle would sediment vertically downward. Despite the fact, that there are 
no direct forces between the particles, they inhuence their motion. This ’’interaction” of the 
immersed particles is mediated by huid and is called the ’’hydrodynamic interaction”. 


gravity field 


r 


r 


Apart from the considerations for the hnite number of particles, Smoluchowski analyzed 
also an inhnite set of particles. He concluded, that behavior of suspension strongly depends 
on the shape of the system. The shape matters, even if its boundaries are extended to 
inhnity. In other words, Smoluchowski identihed the problem of long-range hydrodynamic 
interactions. Another important feature of the hydrodynamic interactions is their many- 
body character. Motion of three particles cannot be described as a superposition of the 
two-particle characteristics. Similar holds for larger number of particles. In general, many- 
body characteristics are needed in the macroscopic considerations for suspensions. From 
the perspective of transport properties, even the two-body hydrodynamic interactions are 
problematic. Analysis of the two-body problem reveals, that two particles at a small distance 
in incompressible, viscous fluid, strongly ’’interact” hydrodynamically. In order to keep 
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constant velocity of the approaching particles, asymptotically an inhnite force is needed [3]. 

Extension of Einstein’s analysis for more concentrated systems appeared to be difficult, 
because of the long-range hydrodynamic interactions. One of the hrst successful approaches 
was made by Saito [4], who obtained the following formula for the effective viscosity, rjes/r] = 
(1 -|- 30/2) /(I — 0). Saito took the hydrodynamic interactions into consideration partially. 
He also discussed the long-range character of the hydrodynamic interactions and strongly 
emphasized difficulties unsolved at that time [5]. The hrst systematic extension of Einstein’s 
work for more concentrated suspensions was performed by Peterson and Fixman [6]. They 
obtained a virial expansion of the effective viscosity up to the second order, which includes 
the two-body hydrodynamic interactions. It was the hrst approach, in which the transport 
coefficient in the second order was given by absolutely convergent integrals. Despite of this 
success, they did not express the transport coefficients by absolutely convergent integrals 
for higher orders of virial expansion. Therefore, the problem with long-range hydrodynamic 
interactions was still not solved at that time. Solution came with the work of Felderhof, Ford 
and Cohen in 1982 [7]. The above authors considered a dielectric system, but their analysis 
can be directly carried over to the physics of suspensions. They proved, that the dielectric 
constant is a local quantity, which does not depend on the shape of the system. Their idea 
is related to the Brown’s approach, who obtained similar result limited to the lowest terms 
in the single-particle polarizability expansion of a dielectric constant [8]. Felderhof, Ford 
and Cohen also gave the microscopic explanation of the famous Clausius-Mossotti formula 
[9], which is an analog of the Saito formula [4] for the effective viscosity in the physics 
of suspensions. It is worth mentioning here the effective medium approaches [10, 11] and 
their extensions including the two-body hydrodynamic interactions in a more accurate way 
[12, 13]. 

Nowadays, the most prominent statistical physics approach to the short-time transport 
properties of suspensions is the Beenakker-Mazur method [14-16], which was developed and 
applied for different suspensions [17-22]. The method gives reasonable results for a wide 
range of volume fractions, but it does not take the two-body hydrodynamic interactions 
fully into account. It is known from virial expansion [23-25] and from numerical simula¬ 
tions, [26, 27] that the two-body hydrodynamic interactions of close particles are essential 
to grasp the dynamics of the system. Therefore, there is still an open problem in the physics 
of suspensions: formulation of a systematic method, which would take the two-body hydro- 
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dynamic interactions fully into consideration and which would give reasonable results for 
at least the intermediate volume fractions, say 0 ~ 25%. Systematic consideration of the 
two-body hydrodynamic interactions in the Beenakker-Mazur method is difficult, because 
the method relies on the expansion of the transport properties in the series of the so-called 
renormalized fluctuations. This series expansion is then truncated in the second order in the 
fluctuations. To consider the full two-body hydrodynamic interactions in the Beenakker- 
Mazur expansion, one needs summation of all orders in the series, which is impossible in 
practice. 

In this article, we develop the approach of Felderhof, Ford and Cohen. As mentioned 
above, they introduced the cluster expansion of the macroscopic characteristics of dispersive 
media such as e.g. the polarizability of dielectric and the effective viscosity of suspension [7]. 
Felderhof, Ford and Cohen also gave the microscopic explanation of the Clausius-Mossotti 
formula for dielectrics (related to the Saito formula in case of suspensions) [9]. Their cluster 
expansion is a starting point of this article. We perform a rigorous analysis of the cluster 
expansion leading to a formula, which we call ring expansion of the macroscopic charac¬ 
teristics. We also generalize the Clausius-Mossotti approximation, basing on a similarity 
between the Felderhof, Ford and Cohen’s cluster expansion and the ring expansion intro¬ 
duced in this article. Using the generalized Clausius-Mossotti approximation, we calculate 
the effective viscosity and the hydrodynamic function (with the translational short-time self¬ 
diffusion and the collective diffusion coefficient) for suspension of monodisperse hard-spheres 
in equilibrium. 

The generalization (renormalization) of the Clausius-Mossotti approximation based on 
the ring expansion introduced in this article is motivated by the results of the virial expansion 
for the effective viscosity and the sedimentation coefficients [24, 25]. One of the dominant 
contributions to the virial expansion for the sedimentation coefficient on the three-body level 
comes from the terms with a virtual overlap of spheres. The idea of resummation of the 
above terms with a virtual overlap of spheres for more dense suspensions was presented to 
the author of this article by Prof. Bogdan Cichocki, to whom the author is very grateful. 

The outline of this article is as follows. In second section, we describe the suspension on 
the microscopic level and discuss the macroscopic characteristics of suspensions. In third 
section, we repeat the analysis of Felderhof, Ford and Cohen leading to the cluster expansion 
of the macroscopic characteristics. In fourth section, we introduce the novel ring expansion 
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of the macroscopic characteristics, which is a rigorous result. The ring expansion is further 
used in hfth section of the article, to introduce a generalization of the Clausius-Mossotti 
approximation. Here, we also present the short-time transport properties calculated by this 
novel method. The generalized Clausius-Mossotti approximation is discussed and its results 
are compared with the results of the numerical simulations and with the Beenakker-Mazur 
method. 


II. MACROSCOPIC PROPERTIES OF SUSPENSIONS 


We consider suspension of hard spheres of radius a in incompressible Newtonian fluid 
of kinematic viscosity rj. We also assume sufficiently slow motion of the fluid and the 
condition of no slip on the surface of immersed particles. As a result, the fluid is described 
by the stationary Stokes equations with the stick boundary conditions [28]. The stationary 
Stokes equations for the problem of the suspension of N spheres, centered at positions 
X = Ri,..., Rat, freely-moving in ambient flow vq (r), under action of external force density 
held fext (r), can be represented in the following integral form [29, 30] 


U, (X, t) = J dV'Mo (R„ r, r') fext (r') 



In the above equations, the particle velocity held Uj (X; r) is dehned inside the particle, i.e. 
for |r — Rj| < a. For hard spheres it has always the following form 


U,(X;r) = W(X) + f7,(^) X (r-R,), for |r-R,|<a, (2) 


with translational Vj and rotational fi, velocity of the particles. Moreover, f, (X; r) describes 
the force density [31-33] acting on the huid by the surface of the particle number i and is 
dehned by 

fj (X; r) = -a (r) ■ hi (r) 5 (|r - Rj] - a), (3) 
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where a represents the stress tensor in the fluid, n* (r) = (r — R*) / |r — Rj| is a vector 
normal to the surface of the sphere i, whereas 6 {x) stands for the one-dimensional Dirac 
delta function. Gq (r) in equations (1) is the Oseen tensor, 

Go(r) = (l + ff)/(87rr/|r|), (4) 

with f = r/ |r|. The Oseen tensor is a Green function of the Stokes equations [34], hence 
flow of the whole suspension v (r) is given by 

N 

V (r) = Vo (r) ^ / d^r'Go (r - r') -f^ (r'). (5) 

i=i 

The equations (1) are linear both in the ambient flow vq and in the external force density 
fext- Therefore, to describe the response operators Mq, M<, M, and M>, it is sufficient and 
convenient to consider special cases of a single particle problem. Mq in the equations (1), 
in the case of the single particle problem in the external force density field fext (r) , and in 
absence of the ambient flow, vq = 0, yields the velocity field of the particle, 

Ui (Ri; r) = j d^r'Mo (Ri, r, r') fe,,* (r'). (6) 

The single particle operator M< gives the particle velocity field Ui, when the particle is 
placed in the ambient flow vq, 

Ui (Ri; t) = J dV'M< (Ri, r, r') vq (r'). (7) 

M(R, r, r'), called the single particle convective friction kernel, yields the force density 
fi (Ri;r) on the surface of the single particle at the position Ri, when it is placed in the 
ambient flow vq (r), 

fi (Ri; r) = j d^r'M {R,, r, r') Vq (r'). (8) 

Finally, M> describes the force density fi (Ri; r) on the surface of the single particle at the 
position Ri, under the action of the external force fext, 



In this article, we investigate the equations (1) mostly without referring to the specific form 
of the response operators Mq, M<, M>, M. For their detail form, we refer the reader to 
the references [30, 35]. 
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To facilitate further analysis of the equations (1), we omit integral variables in those 
equations, writing them in the following form 


U, (X)=Mo (^)fext+M< (^) 


N 


v„ + ^ G„f, (X) 


i=i 


f,(X)=M(^) 


N 


Vo 


Z (X) 


i=i 


M> (i) t,„. 


( 10 ) 


For the position of the particle i in the single particle response operators Mq, M<, M, M>, 
we also use the following abbreviation: f = Rj. Finally, we write the above equations as 
follows [30], 


1 


/ 

p -] 

N 

r 1 

\ 

u, (X) 

= M(i) 


fext 

Uj (Jf) 




V 

Vo 

3 = ^ 


/ 







(11) 


introducing 6x6 dimensional matrices M and G defined by the below equations. 


and 


M (R, r, r') 


Mo(R,r,r') M<(R,r,r') 
M>(R,r,r') M (R, r, r) 


G 


0 0 

0 Go(r,r') 


( 12 ) 


(13) 


A. Scattering series 


To solve the equations (11), several methods can be used. One of the possible approaches 
is the method of reflections [36]. It relies on taking successive iterations of the equation (11) 
which leads to the following formula: 


u. (A') 

= M(i) 

fext 

N 

+ 5^M(!) GM(j) 

fext 

f.(A) 


Vo 

i=i, 

Vo 







N N 


GM(j)GM(i) 


i=i,fc=i, 
jjti k^j 


fext 

Vo 


+ 


(14) 
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The above expression of the force densities, f* (X), and the velocities of the particles, U* (X), 
has a form of a multiple scattering series. It means, that fj (X) and Uj (X) are given by the 
sum of the scattering sequences, for example: 

M (1) GM (2) (15) 


and 

M (1) GM (2) GM (3) GM (2). 


( 16 ) 


As we interpret - each scattering sequence is a superposition of the single-particle scattering 
operators M (i), which ’’scatter” the flow and of Green functions G, which ’’propagates” 
the flow. 

It is convenient and useful to represent the scattering sequences graphically [29]. The 
above two sequences can be represented respectively by 


2 ---p- 

, 2 

1 -d--- 


fo-- 


( 17 ) 


1 -d 


In general, to represent a scattering sequence graphically, we draw horizontal dashed lines 

-. Each line corresponds to a particle in the scattering sequence. Then, reading the 

sequence from left to right, we put successively: the circle O on the dashed line i for the 
operator M (i) and the vertical line | connecting the dashed lines i and j for the Oseen 
tensor G, when it appears in the configuration M [i) GM [j). 

The scattering series plays a major role in our considerations. We denote the scattering 
series by T^ (X): 


T,, (X) = M (^) p + M (^) GM (j) (1 - p) 

N 

+ M{i)GM{k)GM{j) +.... (18) 

k=l 

Therefore, the velocity U* (X) and the force density f* (X) in the expression (14) are given 
by the formula 


U, (X) 


N 

E ^6 (V) 

i=i 


fext 

Vo 


(19) 





B. Macroscopic response 


To describe properties of suspension on the macroscopic level, we consider an ensemble of 
configurations of particles X = Ri,..., Rat, which is described by a probability distribution 
function p{X). We also introduce the average force density dehned by the equation 

/ N \ 


(f(R,r)) 




R 


h; 


(^, r) ) 


( 20 ) 


, 2=1 


and the average particle velocity held 


(U(R,r)) 


N 

E 


i(R-i)Ui(V, r; 


( 21 ) 


\i=l / 

where the three-dimensional Dirac delta function 6 (R — i) = S (R — R*) and the average 
over the probability distribution ([...]) = J d?Ri... f d?R]^ p (X) [...] are used. Averages 
of the equations (19) - multiplied by proper Dirac delta functions - lead to the following 
expression for the average velocity and the average force density, 


(U(R,r)) 

(f(R,r)) 


/ 


d^R'd^r'T (R, r; R', r') 



( 22 ) 


where the averaged scattering series is denoted by T (R, r; R', r') and dehned with the 
formula 

INN \ 

T (R, r; R', r') = ( 5 ] ^ i (R - ») T„ (A', r, rV (R' - j) ) . ( 23 ) 

\i=l i = l / 

Notice, that in the above operator T (R, r; R', r') , the Dirac delta functions hx positions 
of the hrst (i) and the last (j) particle in the scattering series Tjj at the positions R and 
R' respectively. The average how of the suspension (v (r)) is a combination of the ambient 
how vq (r) , in which the particles are immersed and of how generated by the presence of 
the particles 

(v (r)) = Vo (r) j d^R j dVCo (r, r') (f (R, r')), (24) 

which is obtained by averaging the formula (5). We eliminate the how vq from the equations 
(24) and (22), which leads to the formula 


(U) 


fext 

- TG 

(U) 

= T 


. . 


_(v)_ 


. . 


( 25 ) 
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in which we also facilitate the notation by omitting the integral variables. Its subsequent 
iterations lead to a relation of the average particle velocity (U) and the force density (f) to 
the external force density fext and the average flow of suspension (v), 


(U (R, r)> 

_ {f(R.r)> 

which defines operator given by 

= T (1 + GT)~^. 


= / (R, r; R', r' 


fext (r') 

(v (r')) 


(26) 


(27) 


The equation (26) is directly related to the macroscopic properties of the suspension. For 
example, the effective viscosity r]es can be inferred from the relation between the average 
force density (f (R, r)) and the average suspension flow (v (r)) , when no external forces act 
on the particles, fext = 0. The relation between (f (R, r)) and (v (r)) in this situation results 
from the equation (26), after projection it into the lower half of the double vectors [(U), (f)] 
and [fext, (v)]. In order to do that, we introduce a projector Pl defined by 

= (f> , (28) 

with its transposition Pj. After projection, the equation (26) reads 

(f (R, r)) = y d^P'dV'PiT'" (R, r; R', r') Pf (v (r')) . (29) 

If the operator is known, by calculating the following four rank Cartesian tensor 

(R, R') = (r - R)a (R, P R', r') Pf ] (r'-R'), , (30) 

and by symmetrizing it over the first and the second pair of the Cartesian indexes 

-^0^57 ^ i^a/3S'Y (R-, R-O + (R-, R-O + ^a/S'yS (R, R^ + (R, R^) ) (31) 

we obtain the effective viscosity rjes, using the formula [28, 29] 

ta = r, + lim ^ ^ f d^R f d^R'XSir,^ (R, R') , (32) 

Thermodynamic limit limoo is performed in the above equation. 
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Apart from the short-time effective viscosity T^efr, we also consider the short-time wave 
dependent sedimentation coefficient H (q). The sedimentation coefficient describes response 
of the suspension to the external force of the plane wave form, 

Fext (R) = -FoqReexp (-fqR). (33) 


We show in the appendix A, that under the action of the above force, the average transla¬ 
tional velocity of the particles defined by 

N 


(V(R)) = (5]5(R-R,)W (X)), 


(34) 


, 2=1 


in an isotropic and homogeneous suspension, has also a plane wave form. 


(V (R)) = F (g) qReexp (-fqR). 


(35) 


Linearity of the Stokes equations implies, that the coefficient V (g) in the above formula is 
proportional to the force Fq, 

V{q)=H{q) poFo. 


This formula dehnes the wave dependent sedimentation coefficient F[ (g), which is also called 
the hydrodynamic function. The factor po = 1 /( 671770 ) denotes the Stokes coefficient. R (g) 
is a dimensionless function with the property R (g) —)■ 1 in the limit of a diluted suspension, 
i.e. when the volume fraction 0 —)■ 0. As we also discuss in the appendix A, the microscopic 
expression for the hydrodynamic function R (g) has the following form [29] 


H (g) = — ^Tr [ d?R (R) 

ho 3 [g 

where 3x3 matrix Y (R) is dehned by the following equation 


(36) 


F(R- R') = 


(|7ra3)^ 


lim / Rr / dVPuT" (R, r; R', r') R 


' tdT 

u- 


(37) 


The projector Pjj projects on the upper half of the double vectors [(U), (f)], 


pU 


(U) 

(f) 


(U). 


P^ denotes transposition of Pu- 


(38) 
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It is worth noting, that the hydrodynamic function for the zero wave vector, g = 0, 
describes the sedimentation rate, K, of the suspension in a gravity held, 

K = H{q = Q) (39) 

and is also related to the short-time collective diffusion coefficient Dc, 

D, = DoH{q = 0), (40) 

whereas for inhnite wave vector length H [q ^ oo) is related to the short-time self-diffusion 
coefficient Dg, 

Ds = DqH {q ^ oo). (41) 

In both expressions Dq = ksT/ (pnrja) is the diffusion coefficient of a single particle. 

Both, the effective viscosity, and the hydrodynamic function can be inferred from the 
operator. It is shown by the expressions (30-32) for the effective viscosity r^efr and by the 
equations (36-37) for the wave dependent sedimentation coefficient H (q). Therefore, 
becomes the quantity of the main interest in this article. 

III. FELDERHOF, FORD AND COHEN ANALYSIS OF 

In the hrst stage of our analysis of dehned by the equation (26), we follow the idea 
of Felderhof, Ford and Cohen. They obtained the microscopic expression for for the 
dielectric system in the form of a cluster expansion [7]. The application of their idea to the 
physics of suspensions is straightforward, because the governing equations are similar for 
suspensions and dielectrics [29, 37]. To perform the cluster expansion of the operator T^^’^’or 
the basis of the expression (27), Felderhof, Ford and Cohen introduced the cluster expansion 
of the operator T. 

A. Cluster expansion of T 

In the expression (23), Tjj(r,r';X) includes inhnitely many scattering sequences, as 
shown in the formula (18). There are scattering sequences with different number of particles: 
single particle scattering sequences, e.g. 

1-0- ,4 -O- , (42) 
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two-particle scattering sequences, e.g. 


2 -1 


-3- 4--: 

A 

1 C 

'' 

A, 

’ 2 c 

-J- 


( 43 ) 


and scattering sequences with higher number of particles, up to N. The scattering se¬ 
quences with the same number of particles may include different particles. It is noticeable 
in the examples above, where the hrst scattering sequence is between the particles from the 
group C = {1,2}. The second scattering sequence is between the particles from the group 
C = (2,4}. All the scattering sequences (R- — *) Tij (r, r'; X) S (R' — j) can be 

divided, regarding which particles appear in a scattering sequence. To perform this division, 
from all scattering sequences (R — ■*) Tjj (r, r'; X) 5 (R' — j ), we extract only 

the scattering sequences between the particles from the group C: 


TW (R.riR',r'||C) 

N N 

= all s — particle scattering sequences from 

i=l j=l 

which include all particles from s — particle group C. (44) 


The above dehnition allows to represent the cluster expansion of the scattering series as 
follows 

N N N 

J]^i(R-»)T„(r,r'iX)'5(R'-j) = E Y. T'*' (R, r; R, r'| |C). (45) 

i=l i=l S=1 C(ZX,\C\=S 


In the above expression, \C\ stands for the number of particles in the group C, whereas 
'^ccx |c|=s denotes summation over the s-particle groups of particles among X = (1,..., N}. 
Number of such s—particle groups is given by the Newton symbol ('^). 

Average of the equation (45) over the probability distribution function, leads to the cluster 
expansion for the average scattering series T given be the equation (23), 

T(R,r;R',r') = /5^ (R, r; R', r'| |C) V (46) 


, S=1 C(lX,\C\=s 


Since all particles are identical, i.e. the probability distribution p is symmetric with respect 
to interchange of the positions R*, all terms with the same number of particles s in the 
above expression give the same contribution. Therefore, we simplify the last expression, by 
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taking one s—particle group C = {1,..., s} and multiplying it by the factor It yields 


N 


, 5 = 1 


N 


T (R, r; R', r') = ( (R, r; R', r'111... s) 


(47) 


Introducing s-particle distribution functions defined by 

N\ 


n{l.. .s) = 


d'^Rs+i... (TRnP (1... iV), 


(iV-s)! 

we obtain the cluster expansion of the response operator T in the following form, 

N , 


(48) 


T (R. r; R', r') = -T<”> (R, r; R', r'), 


(49) 


5=1 


T'’> (R, r; R', r') = j ... j d^R,n (1 ... s) T<’> (R r; R', r'| 


(50) 


Its thermodynamic limit is achieved by extending of the summation up to iV = cx) and 
performing the thermodynamic limit of the s—particle distribution functions n. From now 
on, we will consider the suspension in the thermodynamic limit. 


B. Nodal line 

To perform the cluster expansion of operator, Felderhof, Ford and Cohen used the 
relation (27), which may be represented in the following form 

Tirr rpQrp ^ xGTGT - . . . . (51) 

Let us look at the second term, i.e. TGT. Representing the T by the cluster expansion 
(49), produces many terms, each of the form 

(52) 

In the expression the scattering sequences between si particles appearing in 

are ’’connected” by the Green function G with the scattering sequences consisted of 
S 2 particles appearing in Altogether, forms si + S 2 -particle scattering 

sequences. The scattering sequences built from the are of a special type, i.e. 

there is a line G connecting a particle from to a particle from Tb 2 ). This line G is the 
only ’’connection” between the particles from T^b and Tb 2 ). It is critical to distinguish the 
lines G, which are the only connections between some groups of the particles in a scattering 
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sequence. Those G are called the nodal lines [29]. It is described by the following examples. 
In the scattering sequence given by expression (15), there is one propagator G. It is a nodal 
line, because it is the only connection between the groups of particles {1} and {2}. In the 
scattering sequence (16), there are three propagators G. The underlined propagator, 

M(1)GM(2)GM(3)GM(2), (53) 


is a nodal line, because it is the only connection between the group {1} and the group {23}. 
In diagrammatic language, the last scattering sequence reads 


3 

2 


n- 




nodal line 


Id— 


(54) 


It is easy to identify a nodal line in diagrammatic language: if cutting a line of a propagator 
G, divides the diagram into two separate pieces, then the propagator G is a nodal line. 


C. Cluster expansion of T with nodal lines specified 

In the previous section, we indicated, that an important element of the analysis of is 
the notion of the nodal line. Therefore, we perform further division of the scattering series 
T, by specifying the nodal lines in the scattering sequences. 

In agreement with the dehnition (44), T*^®^ (R, r; R', r'| |1... s) represents inhnitely many 
s—particle scattering sequences. We divide them into disjoint sets, by specifying the number 
of the nodal lines and by specifying the groups of particles separated by the nodal lines in the 
scattering sequence. Those sets are characterized by the number of groups g, the number of 
particles in each group Si = |Gi| ,... ,Sg = \Cg\, and by saying which particles from (1 ,..., s} 
are in the group G*. The number of the groups g is larger by one than the number of the 
nodal lines. Since we consider the s-particle scattering sequences, we have the condition si + 
... + Sg = s. We extract from all s—particle scattering sequences T*^®) (R, r; R", r"| |1. .. s) , 
the scattering sequences with specihed groups of particles Gi,...,Gg separated by nodal 
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lines, defining 


T'*l(R,riR',r'||Ci|...|C,) 

= all scattering sequences from (R, r; R', r'| |Ci... Cg) with g — 1 nodal lines 
separating particles from the groups Ci,... ,Cg. (55) 

With the above dehnition, the s—particle scattering sequences can be divided as follows, 
T<*>(R,riR',r'||C) = ^ T'*'(R, r; R', r'||Ci|... |C,). (56) 

g=l si+...+Sg=s Ci,...,Cg(ZC, 

\Cl\ + ...+ \Cg\=S 

Here, Ci,...,CgCC, denotes summation over all possible divisions of the set of s particles 
s} into g subsets, with si particles in the first subset, S 2 particles in the second 
subset, etc. There are s\/ (si!.. .s^!) possible divisions. 

Let us now consider the lowest order term (C) in the expression (56), i.e. the term 
with (7 = 1 , which have no nodal lines. The dehnition (55) implies that 

(R, r; R', r'| |C) = all scattering sequences from (R, r; R', r'| |C) without nodal lines. 

The above scattering sequences without nodal lines play signihcant role. They are called 
the irreducible scattering sequences [29]. 

The second order term in the expression (56) is the term with g = 2, 

Y, Y1 T<*>(R,riR',r'||Ci|C 2 ). (57) 

Si-\-S2=S Ci,(72CC, 
iCil + lCsN^ 

From the dehnition (55), it follows, that {Ci\C 2 ) has one nodal line G separating the 
particles from the groups Ci and C 2 . Therefore, all scattering sequences in (G 1 IG 2 ) 
have the following structure: hrst - there are some rehections between the particles from the 
group Cl, then - there is exactly one rehection G to a particle in the group C 2 (nodal line), 
and then - there are rehections between the particles from the group C 2 . The rehections 
between the particles from the group Ci must be irreducible (without a nodal line). The 
same holds for the group G 2 . It suggests, that (G 1 IG 2 ) has the following structure, 

(GilGs) = {Cl) GT^'"^ {C 2 ), (58) 


16 



in which the irreducible scattering sequences {C) appear and the nodal line is written 
explicitly. The above formula can be simply proved using the dehnition (55). Similar results 
hold for the higher terms of the expansion (56), 

tW (Cil... \Cg) = (Cl) G ... {Cg ), (59) 

for the groups Ci,..., Cg, including si,..., Sg particles respectively. In the above formula, 
the nodal lines separate different irreducible sections (Q) of the scattering sequences 
(Ci| ... \Cg). Each irreducible section T^®^ (Q) is referred to as ’’the block”. Moreover, 
by ’’the block structure” - we mean the way particles are distributed in the blocks. A block 
structure is specihed as follows: Ci|... \Cg. The block structure of the scattering sequence 
given by the expression (16) is 1|23. It is convenient to introduce the following notation for 
the irreducible scattering sequences, 

S, (R, r; R', r'| |C) = 1 t<‘I (R. r; R', r'| |C), (60) 

with the factor s! Using the above two formulas, we rewrite the expansion (56) as follows 

T<*)(C) = 5; ^ ^ s,!...s,!S,(C,)G,,,GSaC,), (61) 

|Cl| + ... + |Cg|=S 

The above representation used in the cluster expansion of T operator represented by equa¬ 
tions (49) and (50), after a simple algebra, leads to the expression 

oo „ 

T = E E / dCi... dCg n(Gi,..., G3 )S7 (Gi)G ... GS/(GJ. (62) 

g=l Ci,...,C 9 

Cancellation of the factors sj results from the symmetry of the probability distribution 
density p and the fact, that the particles are identical. The symbol Yhci Cg I • • • / 
denotes summation over different numbers of particles in each of g groups and integration 
over the positions of particles as follows 

f dC,...dCgfiC^,...,Cg) = 

OO oo « 

J2 T, /(R!,...,R;,,...,Rf,...,RS,) (63) 

ni = l ^5=1 ^ 
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D. Cluster expansion of 

We go back to the expression (51), 

Tirr rpQrp ^ TGTGT - ..., (64) 

in order to derive the cluster expansion of the operator. As we discussed before, the 
average scattering series T includes all possible scattering sequences. TGT in the equation 
(64) produces scattering sequences with at least one nodal line, TGTGT with at least two 
nodal lines, etc. The analysis of the above equation relies on a consideration of the scattering 
sequences with given number of nodal lines. Therefore, we will consider terms with different 
block structures Gi| ... |Gg on the right-hand side of the equation (64). 

Let us focus hrst on the block structure Ci, i.e. on the scattering sequences without 
nodal lines. They appear only in the first term, T, of the equation (64), because the higher 
terms, TGT, TGTGT, ..., include at least one nodal line. In the expression (62) for T, 
the irreducible scattering sequences come only from the term g = 1. Therefore, all terms on 
the right hand-side of the equation (64) with the block structure Ci, are of the form 

n(C'i)S,(C'i). (65) 

Next, we consider the terms on the right-hand side of the equation (64) with the block struc¬ 
ture G 1 IG 2 , i.e. with one nodal line. Such scattering sequences, i.e. S/(Gi)GS/(G 2 ), appear 
only in the hrst term, T, and in the second term, TGT. A contribution from the T comes 
from the term with g = 2 of the equation (62) and is of the form n {C 1 C 2 ) S/(Gi)GS/(G 2 ). 

A contribution from TGT has a form of n (Gi) n (G 2 ) S/(Gi)GS/(G 2 ) and is produced by 
the terms with g = 1 in both T. Those two terms altogether yield 

(n (Gi, G 2 ) - n (Gi) n (G 2 )) S,(Gi)GS,(G 2 ). ( 66 ) 

For the block structure consisted of the three groups G 1 IG 2 IG 3 , analysis of the right-hand 
side of the expression (64) leads to the contribution of the following form 

(n (G 1 G 2 G 3 ) - n (G 1 G 2 ) n (G 3 ) - n (Gi) n (G 2 G 3 ) + n (Gi) n (G 2 ) n (G 3 )) S,(Gi)GS,(G 2 )GS,(G 3 ). 

(67) 
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In such manner, analysis of all block structures Ci \ .. .\Cg is possible. Functions ap¬ 
pearing along with the block structures of the form S/(Ci)G ... GSi{Cg), are denoted by 
b (Gil... \Cg) and are called the block distribution functions [7]. Therefore we have 

6 (Gi)=n(Gi), ( 68 ) 

b (Gi IG 2 ) = n (Gi, G 2 ) - n (Gi) n (G 2 ), (69) 

b (G 1 IG 2 IG 3 ) = n (G 1 G 2 G 3 ) - n (G 1 G 2 ) n (G 3 ) - n (Gi) n (G 2 G 3 ) + n (Gi) n (G 2 ) n (G 3 ), 

(70) 

for the block structures up to three groups. Expressions for the block distribution func¬ 
tions for higher number of groups are more and more complicated. Nevertheless, the block 
distribution functions may be calculated from the following recursive formula [29, 38] 

6 (G) = n(G) (71a) 

6 (Gi|... \CkCk+i \... \Cg) = &(Gi|... \Cu)b{Ck+i \... IG^) + 

+ b{C^\...\Ck\Ck+i\...\Cg). (71b) 

The above analysis leads to the following cluster expansion of operator: 

= E E / dGi... dCgb{C,\... |G,)S,(Gi)G... GSi{Cg). (72) 

g=l Ci,...,C 3 

It is worth noting, that the structure of operator is similar to the structure of T operator 
given by the expression (62). The only difference lies in the distribution functions: in - 
the block distribution functions b (Gi| ... \Cg) appear, whereas in T operator - the standard 
s-particle distribution functions n (Gi... Cg) appear. 

E. Self scattering sequences 

There are phenomena in the physics of suspensions, in which only a part of the scattering 
series T, given by equation (23), plays a role. An example of the above is the self-diffusion 
coefficient. It is related only to those scattering sequences in T, which start and end at the 
same particle. The scattering sequences, which start and end at the same particle, we call the 
self-scattering sequences. The self-scattering sequences are irreducible, because there are no 
nodal lines in any scattering sequence, which starts and ends at the same particle. Therefore, 
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the self-scattering sequences are related to the irreducible scattering sequences S/, as 
follows, 

(R, r; R', r'| |C) = all scattering sequences from S/(R, r; R', r'| |C), 

which start and end at the same particle. (73) 

The average of the self-scattering sequences will be denoted by B, 

B(R,r;R',rO = ^ / dCi n(C'i)Sf ^(R, r; R', r'| ICi). (74) 

Cl 

With the above quantity, the scattering sequences in T can be divided into the self-scattering 
sequences and the scattering sequences, which start and end at different particles (off- 
scattering sequences). The former are given by B (R, r; R', r'), the latter are denoted by 
Toff. Therefore, 

T = B + Toff. (75) 

Operator can be divided in a similar manner, 

= B + TjCi^. (76) 

As shown in the above formulas, the self-part of both operators T and T'^''' are the same. 

It is worth noting, that the self-scattering sequences B can be calculated from off- 
scattering sequences To® as follows [39] 

B (R, r; R', r') = mS (R - R') M (R, r, r') +5 (R - R') [ToffGM] (R, r; R, r). (77) 

The second term in the above expression corresponds to scattering sequences of the following 
structure. The sequences start at the particle at position R, then go with all possible 
scattering sequences Toff to some particle, and then they come back with one reflection G 
to the starting particle at R. Moreover, in the above equation ni stands for the one-particle 
distribution function. 

IV. RING EXPANSION OF 

The cluster expansion (72) of the response operator was derived by Felderhof, Ford 
and Cohen about three decades ago [7]. This achievement allows to express the transport 
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coefficients of dispersive media, such as suspensions and dielectrics, by absolutely convergent 
integrals. Despite this important step done by the three scientists, no reasonable statistical 
physics method of calculation of transport properties of suspensions can be found in current 
literature. We have already discussed this point in the introduction. 

In order to give a motivation of our approach introduced further in this article, we invoke 
the effective Green function Gefr dehned by 

Geff = G + GTG. (78) 

In the case of homogeneous system in the thermodynamic limit, the effective Green function 
Geff of suspension is related to the Green function G of pure fluid according to the following 
formula, [29] 

Geff (r, r') ^ (r — r'), for |r — r'| —)■ cx). (79) 

VeS ( 0 ) 

The above formula describes the asymptotic decay of the effective propagator Geff (r, r') for 
the large distances between the points r and r'. The decay is governed by the effective 
viscosity of suspension ? 7 eff (0), which depends on the volume fraction 0 = 47 rnia^/ 3 , with 
Hi standing for the number density of the particles in suspension. 

Let us consider relation (27) between T and T"’’’, which can be inverted and then repre¬ 
sented in the following way, 

T = + T'^GeffT'", (80) 

where the formula for the effective propagator Geff = G (l — T^'’'’G) ^ is used. Gonsidering 
only the off-scattering sequences in the above expression, we receive the equation 

Toff= T*" + T'''''GeffT '''0 (81) 

obtained with the application of the formulas (76) and (75). Let us notice, that Toff (R, r; R', r') 
0 for overlapping conhgurations, i.e. for |R — R'| < 2a. It results from the assumption, that 
the hard spheres in suspension cannot overlap and in the expression (62), the distribution 
function n (C) vanishes for the overlapping conhgurations. The vanishing of Toff (R, r; R', r') 
for the overlapping conhgurations has a consequence in the last equation. It reduces to 

/ (R - R') T*- (R, r; R', r') = -/ (R - R') [T^GeffT-] (R, r; R', r'), (82) 
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after its multiplication by a function / (R — R') defined by 


, 1 for |R — R'l < 2a 

/(R-R') = <^ 

0 for |R — R'l > 2a 


( 83 ) 


which equals 0 for nonoverlapping and 1 for overlapping conhgurations of the two particles 
centered at R and R'. 

The equation (82) shows, that in there are some contributions with the effective 
propagator Geff. It suggests, that the propagator G appearing in the cluster expansion (72) 
of Felderhof, Ford and Cohen can be renormalized. In other words, T*'’'’ can be given by the 
following formula. 


T 


irr 


E E dC,... dC,H{C,\... |C',)S,(C'i)Geff... GeffSKG,), 

d=l Ci-Cd'' 


(84) 


with yet unknown functions H (Gil ... \Cd ), which we call the block correlation functions. 
The above expression for T'’'’’ has the same structure as expression (72), but contains the 
effective propagator Ges, instead of the propagator G and also contains the block correlation 
functions H (Gi|... \Cd) , instead of the block distribution functions b (Gi| ... \Cg). We call 
the expression (84) for T''’’' operator, the ring expansion, in order to differentiate it from the 
cluster expansion (72) of this operator introduced by Felderhof, Ford and Cohen. 

Below, we prove the ring expansion (84) and derive a formula for the block correlation 
functions H {Ci \.. .\Cd)- We will use the similar approach, as in the derivation of the 
Felderhof, Ford and Cohen’s formula (72) in the previous section. We will consider the 
right-hand sides of both expressions (72) and (84), considering scattering sequences with a 
given block structure Gi|... \Cg. 

We start with the block structure S/(Gi), i.e. the block structure without a nodal line. 
In the expression (84), the block structure S/(Gi) without a nodal line appears only in the 
term d = 1. In the expression (72), the block structure S/(Gi) also appears in the lowest 
order term g = 1 only. Therefore, equality of the expressions (72) and (84), on the level of 
the irreducible scattering sequences S/(Gi), is possible if we assume 


6(Gi) = R(Gi). 


(85) 


Before further considerations for a general block structure, let us hnd all terms in the 
ring expansion (84), which have the block structure G 1 IG 2 IG 3 IG 4 . For the purpose of the 
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above, we need the cluster expansion of the effective Green function Geg, 


dCi... dCg n(C'i,..., C'JGSKC'OG ... GS/(Gg)G, (86) 

s=o Cu-,Cg'^ 

which is a straightforward consequence of the expressions (78) and (62). In the above 
formula, the term with g = 0 corresponds to the Oseen tensor G. The effective Green 
function introduces one, two, and more nodal lines in the block structure S/ (Gj) Ge//S/ {Cj). 
Therefore, the scattering sequences S/ (Gi) GSj (C2) GSj (G3) GS/ (G4) appear in the ring 
expansion (84) in the term corresponding to d = 2, 

//(G,JG,JS,(G,JGeffS,(G,J, (87) 


in the term d = 3, 


dd(G,JG,,|G,3)S,(G,JGeffS,(G,JGeffS,(G,3), 


( 88 ) 


and in the term d = 4, 


/d(G,JG,,|G, 3 |G,JS,(ajGeffS,(G,JGeffS,(G, 3 )GeffS,(aj. (89) 

Other terms in the ring expansion (84) do not contain the scattering sequences with three 
nodal lines, because the term d = 1 contains no nodal lines, whereas the terms corresponding 
to d > 5 contain at least four nodal lines. 

Each of the terms in the expressions (87-89) contains many different scattering sequences 
produced by the cluster expansion ( 86 ) of the effective Green function. In the case of the 
expression (87), only the term n (G2G3) GS/ (G2) GS/ (G3) G in the expansion ( 86 ) produces 
the block structure G 1 IG 2 IG 3 IG 4 , yielding 

dd (G1IG4) n (G2G3) S/ (Gi) GS/ (G2) GS/ (G3) GS/ (G4) . (90) 

In order to obtain this, we assume ii = 1, f 2 = 4 in the expression (87). In the case of 
the expression ( 88 ), there are two possibilities leading to the block structure G 1 IG 2 IG 3 IG 4 . 
The first possibility corresponds to the situation, when the hrst (left) propagator Geg in the 
expression ( 88 ) introduces one nodal line and the second propagator introduces two nodal 
lines. Therefore, we assume ii = 1, *2 = 2, *3 = 4 and the contribution of the term given by 
the expression ( 88 ) is 

dd (G1IG2IG4) n (G3) S/ (Gi) GS/ (G2) GS/ (G3) GS/ (G4) . (91) 
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The second possibility corresponds to an opposite situation, when the first propagator in 
expression ( 88 ) introduces two nodal lines and the second propagator introduces one nodal 
line. Here, we assume ii = 1, i 2 = 3, is = 4 and the contribution is 

H (C'i|C' 3 |C' 4 ) n {C 2 ) S/ (Cl) GS, (^ 2 ) GS, (Cs) GS, (^ 4 ). (92) 

In the case of expression (89), there is only one possibility to obtain the block structure 
C 1 IC 2 IC 3 IC 4 , i.e. when all propagators Ges introduce only one nodal line G. In this case, 
we have ii = 1, *2 = 2, is = 3, *4 = 4 and obtain the following contribution, 

H (C'i|C' 2 |C' 3 |C' 4 ) S, (Cl) GS, (C 2 ) GS, (C 3 ) GS, (C 4 ). (93) 

Finally, all terms in the expression (84), which have the block structure C 1 IC 2 IC 3 IC 4 (i.e. 
terms given by expressions (90-93)), after comparison with the term containing the same 
block structure C 1 IC 2 IC 3 IC 4 from the equation (72), lead to equality 

b (C 1 IC 2 IC 3 IC 4 ) = H (C 1 IC 4 ) n (C 2 C 3 ) + H (C 1 IC 2 IC 4 ) n (C 3 ) 

+ H (C 1 IC 3 IC 4 ) n (C 2 ) + H (C 1 IC 2 IC 3 IC 4 ). (94) 

Similar considerations for the block structures C 1 IC 2 and C 1 IC 2 IC 3 lead to the expressions 

6 (Ci|C 2 ) = 77(Ci|C2), (95) 

b (C 1 IC 2 IC 3 ) = H (C 1 IC 2 IC 3 ) + H (C 1 IC 3 ) n (C 2 ). (96) 

The above considerations for the block structure C 1 IC 2 IC 3 IC 4 , leading to the formula (94), 
can be generalized to the case of a block structure Ci| ... |Cg consisted oi g >2 groups. The 
block structures Ci| ... \Cg appear in the ring expansion (84) in the terms d = 2,..., g only, 
i.e. the terms of the form S/ (Ci) Ge// ... Ge//S/ (C^). As the term ( 88 ) for the case <7 = 4 
introduces the block structure C 1 IC 2 IC 3 IC 4 in two ways, each of S/ (Ci) Ge// ... Ge//S/ (C^) 
can introduce the block structure Ci \ .. .\Cg in several ways. All terms can be uniquely 
classihed by specihcation of d groups QJ ... \Ci^ among Gi| ... \Cg, which come from the 
blocks S/ in the expression S 7 (Gi)Ge// ... Geff^iiCd)- The edge groups must be the same, 
therefore A = 1 and id = g. Each set of numbers 1 = ii < A < • • • < A-i < id = 9 
corresponds to a single term in the expression S/(Gi)Ge// ... Ge//S/(Grf), which produces 
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the block structure S/(Ci)G... GS/(Cg). Comparison of all terms in the expansion (84) 
and (72), producing the block structure S/(Gi)G ... GSi{Cg) yields 

MCi|...|c,) = 

d=2 l=ii<i2<---<id=9 

n({C„ ... C.J \ {C..C.J) ...a ... O,} \ }) , (97) 

which is valid for g > 2. The symbol ’\’ denotes a difference of sets of the particles, e.g. 
{12567} \ {56} = {127}. We assume that for the empty set n(0) = 1. The above formula 
is a recursive expression for the block correlation functions H. 

A. Comparison of ring and cluster expansion 

The ring expansion (84) introduced in the previous section is a rigorous expression for 
the response operator This is an alternative formula to the cluster expansion (72) of 
Felderhof, Ford and Cohen. The cluster expansion and the ring expansion have the same 
structure: the irreducible scattering sequences S/ connected by the propagators (G or Geg), 
are averaged over conhgurations of particles, weighted with the distribution functions {b or 
H). From that perspective, and due to fact, that the effective Green function Gefr appears 
in the ring expansion instead of G, the formula (84) can be seen as the renormalized cluster 
expansion. 

There are two important differences between the ring and the cluster expansion. The hrst 
difference lies in the propagators. The propagator G, which appears in the cluster expansion 
(72), includes only information concerning liquid. On the other hand, the effective Green 
function Gefr in the expression (84), contains macroscopic information about suspension. It 
is exhibited by the appearance of the effective viscosity ?7efr in the asymptotic form of the 
effective propagator for large distances showed in the formula (79). The effective viscosity 
r]es of hard-sphere suspension may signihcantly differ from the viscosity r] of pure liquid - 
especially for the higher volume fractions 0. 

The second difference between the ring and the cluster expansion lies in the distribution 
functions. The block distribution functions b appearing in the cluster expansion are given 
with the formula (71). The block correlation functions H, which appeared in the ring 
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expansion, are defined with the expression (97). There is an essential difference between b 
and H. It is related to the cluster property of the distribution functions n, which we assume 
in this article. The cluster property relies on the factorization of the distribution function 
n (C1C2) in the limit of large distance between the groups Ci and C2, 

^(C'lC's) ^n(C'i)n(C'2). (98) 

A use of the above cluster property of the distribution functions n in the equation (70), 
when the group C2 in the middle of the block structure C1IC2IC3 goes away from the other 
groups, leads to the following factorization of the block distribution function 

h{C,\C2m ^h{C,\C^)h{C^). (99) 

The cluster property of the distribution functions n applied in the equation (96) for 
77 (Cl I(7*21 Cs) in the same limit - when the group C 2 goes away - results in the follow¬ 
ing decay of the block correlation function, 

77(C'i|C'2|C'3) ^ 0. (100) 

The above asymptotic decay is a motivation for the name of the block correlation func¬ 
tions 77. The above property is also a motivation for the name of the expression (84) - 
i.e. ring expansion. Two subsequent blocks, S/ (C*) and S/ (Cj+i), in the ring expansion 
(84) are ’connected’ by, both, the effective propagator Ges and by the correlation function 
77 (... |Cj|Cj+i|...) - both ’connections’ vanish, when C* goes away from Cj+i. We imagine, 
that such double-connection of the S/ {Ci) and S/ (Cj+i) form a ’ring’. 

The ring expansion of represented by the equation (84) along with the expression 
(97) for the block correlation functions is the main analitycal result of this article. It is an 
alternative expression to the Felderhof, Ford and Cohen’s cluster expansion represented by 
the formula (72). Our ring expansion appears as a result of a resummation performed on the 
level of the cluster expansion. This resummation procedure leads to the ring expansion which 
has similar structure as the structure of the cluster expansion. The role of the Oseen tensors 
in the cluster expansion - after resummation - is played be the effective Green function Gefr 
given by the formula (78). The effective Green function has a physical interpretation because 
it relates the force (generating the ambient flow) with the velocity held of the suspension - in 
contrast to the Oseen tensor which relates the force with the velocity held of a pure liquid. 
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Moreover, the effective Green function is related to the effective visicosity, as the expression 
(79) shows. Because the resummation procedure leads to a similar structure as the structure 
of the starting expression, we call this procedure the renormalization. Consequently, the 
effective Green function may also be called the renormalized (effective) Green function. The 
ring expansion is further used in the next section to introduce a new method of calculations 
of transport properties of suspensions. 

The above derivation of the ring expansion is presented in the short-hand notation which 
emphasizes the idea of the underlying physics. It is worth presenting the ring expansion 
without the short-hand notation. Following the expression (63) the lowest two terms of the 
ring expansion are given by 


T'’-'- (R, r; R', r') = 

OO „ 

ni=l 

00 00 « 

EE/ d^R \... ... d^Rl^d\''dh-'''d^R''(fR'''H (R},..., R^JR^,..., R^ J 

rii=l 712 = 1 


X S,(R, r; R", r"| IRJ,..., R'JGeff (r", r'") S,(R'", r'"; R', r'| |R2, ..., Ri 


) 

712 / 


( 101 ) 


V. RENORMALIZATION OF CLAUSIUS-MOSSOTTI APPROXIMATION 

After the derivation of the cluster expansion (72) of the operator, Felderhof, Ford 
and Cohen gave the microscopic explanation of the Clausius-Mossotti formula [9]. It is an 
expression for the relative dielectric constant of a nonpolar dielectric system. It may be 
derived using the macroscopic considerations [40]. Felderhof, Ford and Cohen explained the 
Clausius-Mossotti formula on the microscopic level, showing a class of terms in the cluster 
expansion of T'*''', which leads to the Clausius-Mossotti relation. 

Going along the line of the explanations of Felderhof, Ford, and Cohen [9] for a dielectric 
system, the Clausius-Mossotti relation is obtained, when the operator dehned with the 
following formula 

(1 + [hG] T'’-'-) , (102) 

is approximated by the single particle term, 

(R, r; R, r') ^ md (R - R') M (R, r, r'). (103) 
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In the above definition (102), a superposition between the quantities (R, r; R', r') and 
[hG] (R, r; R', r') = h (R, R') G (r, r') appears, h (R, R') stands for the two-particle corre¬ 
lation function. It is worth noting, that in the reference [9], instead of the two-particle cor¬ 
relation function h (R, R'), function —/ (R, R') with / defined by the formula (83) appears. 
Both possibilities lead to the Clausius-Mossotti relation. We call the operator TgJvi - the 
Clausius-Mossotti operator - because it is a straightforward generalization of the Clausius- 
Mossotti function. In case of suspensions, the above procedure leads to the Saito formula 
for the effective viscosity [41], namely 


VeS _ 1 + |0 
7 ] 1 — 0 


(104) 


The Clausius-Mossotti approximation is expressed by the approximated formula (103), 
applied to the Clausius-Mossotti operator, defined by the equation (102). In this equation, 
the Oseen tensor G appears. In the previous section, we show that the propagators G in the 
cluster expansion (72) can be renormalized and, as a result, operator can be represented 
by the ring expansion (84), with the effective Green function Ges appearing instead of G. 
This suggests to define renormalized Clausius-Mossotti operator as follows. 


Ti;fcM = T"(l + |AG„lT") 


(105) 


in analogy to the formula (102), and to generalize the Clausius-Mossotti approximation by 


-*-RCM ~ 


B, 


(106) 


in analogy to the approximation (103). In the latter equation, instead of the single parti¬ 
cle response operator M appearing in the Clausius-Mossotti approximation, we take into 
account all self-scattering sequences B, introduced before with the formula (74). The equa¬ 
tions (105), (27), (78), (75), (77) along with the approximation (106) define the renormalized 
Clausius-Mossotti approximation. Those equations form a close system of equations for op¬ 
erators B, Geff, T, and Tog. The system can be solved for given volume fraction 

0 (or the single particle density ni) and for given two-particle correlation function h (R, R'). 
We solve those equations numerically. Not to interrupt our line of reasoning, we refer the 
reader to the appendix B containing the technical details of our numerical calculations. From 
found within the renormalized Clausius-Mossotti approximation, one can calculate fur¬ 
ther the short-time transport characteristics, such as the effective viscosity ? 7 efr from the 



equations (30-32) and the wave dependent hydrodynamic function H (g) with the collective 
diffusion Dc and the self-diffusion Dg coefficient from the equations (36-41). 

Before presenting in the next section the results for the transport characteristics calcu¬ 
lated within the generalized Clausius-Mossotti approximation, in what follows we comment 
on its physical meaning. To this end we discuss the T operator (see equation (22)) obtained 
within the Clausius-Mossotti approximation denoted by . The Clausius-Mossotti ap¬ 
proximation (103) - by equations (102) and (27) - leads to the expression 

T^^ = niM[l - [c/G]niM]"\ (107) 

where g is the radial distribution function related to the correlation function h as follows, 
g = 1 + h. The above formula for T can be interpreted in terms of scattering sequences. 
The T operator - which on rigorous level is given by the sum of all possible scattering 
sequences as the equation (18) shows - in the Clausius-Mossotti approximation is given by 
a sum of scattering sequences in which the reflections never go back to a particle which 
already reflected the flow. Moreover, there are correlations g only between neighboring 
particles in the scattering sequences. On the other hand the generalized Clausius-Mossotti 
approximation (106) - along with the rigorous equations (105), (78) and (27) - leads to the 
following formula for the T operator (denoted by 

rj^RCM = B [1 - {[gG] + h[GTG]) B]"^. (108) 

The above equation differs from the expression (107), in particular by the term h[GTG]. 
This term contains the dominant terms in the virial expansion on three-particle level for 
the sedimemntation coefficient (cf. 64 coefficient in the reference [24]) and for the effective 
viscosity (cf. z /2 coefficient in the reference [25]). Therefore we expect that the renormal¬ 
ized (generalized) Clausius-Mossotti approximation will give more accurate results than the 
original Clausius-Mossotti approximation. 

VI. RESULTS AND DISCUSSION 

The renormalized Clausius-Mossotti approximation introduced in the previous paragraph 
allows to calculate the short-time transport properties of suspension, when the volume 
fraction 0 and the two-body correlation function h (R, R') are given. Within the renor¬ 
malized Clausius-Mossotti approximation, we perform calculations for the volume fractions 
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FIG. 1. Inverse of the translational short-time self-diffusion coefficient Dg, equation (41), normal¬ 
ized by the self-diffusion coefficient of a single particle, Dq = ksT/ (Gvrr/a), as a function of volume 
fraction cj) for monodisperse suspension of hard-spheres in equilibrium. Black (solid line) - the 
renormalized Clausius-Mossotti approximation introduced in this article, red (long-dashed) line - 
numerical simulations of Abade at al. [43], blue (short-dashed) line - the revised Beenakker-Mazur 
method [44]. 

0 = 0.05,0.15,0.25,0.35, and 0.45. For each volume fraction, we use the two-particle corre¬ 
lation function in the Percus-Yevick approximation for the hard-sphere potential [42]. Our 
results are presented in hgs. 1- 5, which show the translational short-time self-diffusion 
coefficient Dg^ the sedimentation coefficient K, the hydrodynamic function H{q), and the 
effective viscosity coefficient rjes, respectively. 

In the figs. 1- 5 we also present results of the numerical simulations and the (revised) 
second order Beenakker-Mazur method [44], which is nowadays the most comprehensive 
theoretical scheme for calculations of the short-time transport properties of suspensions. At 
this point it is worth noting, that it is difficult to realize experimentally a monodisperse 
suspension of hard-spheres [46] and to measure its volume fraction precisely [47]. However, 
if the experimental conditions satisfy the assumptions underlying our theoretical model 
(such as monodisperse, hard-sphere potential, regime of the zero Reynolds number), then 
this suspension would have characteristics consistent with the precise numerical simulations. 
Therefore, in this article, we assess the theoretical methods for monodisperse hard-sphere 
suspensions by comparing them with numerical simulations, instead of comparing with the 
experimental works [48-56]. 
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FIG. 2. Inverse of the sedimentation coefficient equation (39), as a function of volume fraction 
(j) for monodisperse suspension of hard-spheres in equilibrium. Black (solid line) - the renormalized 
Clausius-Mossotti approximation introduced in this article, red (long-dashed) line - numerical 
simulations of Abade at al. [43], blue (short-dashed) line - the revised Beenakker-Mazur method 
[44], green (dot-dashed) line - Clausius-Mossotti approximation defined by the equation (103). 



FIG. 3. The hydrodynamic function H (q), equation (36), as a function of wave vector for volume 
fractions (j) = 0.05, (j) = 0.15 and (p = 0.25, for monodisperse suspension of hard-spheres in equi¬ 
librium. Black (solid line) - the renormalized Clausius-Mossotti approximation introduced in this 
article, red (long-dashed) line - numerical simulations [43], blue (short-dashed) line - the revised 
Beenakker-Mazur method [44]. 

Our results for the effective viscosity and the hydrodynamic function H (q) (with its 
low and high q behavior) presented in hgs. 1-5 obtained within the renormalized Clausius- 
Mossotti approximation, when compared with the numerical simulations and the revised 
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FIG. 4. The hydrodynamic function H {q), equation (36), as a function of wave vector for volume 
fractions cj) = 0.35 and (f) = 0.45, for monodisperse suspension of hard-spheres in equilibrium. Black 
(solid line) - the renormalized Clausius-Mossotti approximation introduced in this article, red (long- 
dashed) line - numerical simulations [43], blue (short-dashed) line - the revised Beenakker-Mazur 
method [44], 



FIG. 5. The relative effective viscosity ?7efr/?/ (high frequency, low shear), equation (32), as a 
function of volume fraction (p for monodisperse suspension of hard-spheres in equilibrium. Black 
(solid line) - the renormalized Clausius-Mossotti approximation introduced in this article, red 
(long-dashed) line - numerical simulations of Ladd [45], blue (short-dashed) line - the revised 
Beenakker-Mazur method [44], green (dot-dashed) line - Saito formula [4]. 

second order Beenakker-Mazur method - can be summarized as follows. For the volume 
fractions 0 < 30%, the relative error (with respect to the numerical simulations) of the renor¬ 
malized Clausius-Mossotii approximation is a few times less or comparable with the relative 
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error of the Beenakker-Mazur method - it holds for the effective viscosity, the sedimentation 
coefficient, the self-diffusion coefficient, and for almost whole range of the wave vectors q of 
the hydrodynamic function. The only exception is a range of the wave vectors 2 < 2aq < 5 
placed between g = 0 and q in the principal maximum of the hydrodynamic function H (g) 
(vide hgs. 3 and 4). Here, the Beenakker-Mazur method is consistent with numerical simu¬ 
lations. For volume fractions above 0 ~ 30%, the Beenakker-Mazur method leads to better 
agreement with numerical simulations than the renormalized Clausius-Mossotti approxima¬ 
tion, for the effective viscosity and for most of the wave vectors g of hydrodynamic function. 

Before comparison of the results of our renormalized Clausius-Mossotti approximation 
dehned by the equation (106) with the Clausius-Mossotti approximation dehned by the 
formula (103), it has to be emphasized that the term Clausius-Mossotti approximation may 
refer to two different variants of approximation. First variant is given by the equations (103) 
and (102), and is used in this article. The second variant is also given by the equations 
(103) and (102), but instead of the two-body correlation function h, its lowest virial term, 
i.e. the Mayer function for hard spheres [57], is used. The Clausius-Mossotti approximation 
in case of the effective viscosity leads to the Saito formula (104), whereas in case of the 
hydrodynamic function denoted in this approximation by HcM{q) - because of the fact that 
for the large wave vectors hydrodynamic function is related to the self-scattering sequences 
and the Clausius-Mossotti approximation takes into account only one single particle-term 
among self-scattering sequences - it gives the following result 

lim HcM{q) = 1- (109) 

g—>-oo 

Therefore the self-diffusion coefficient in the Clausius-Mossotti approximation does not de¬ 
pend on the volume fraction of suspension. The opposite limit of the zero wave vector related 
to the sedimentation coefficient is presented in the £g. 2. 

It is worth shedding light on the results of our article - the derivation of the ring expan¬ 
sion (84) and the formulation of the renormalized Clausius-Mossotti approximation - from 
the perspective of the hydrodynamic interactions and statistical physics. The Beenakker- 
Mazur method is currently the most comprehensive statistical physics method to calculate 
the short-time transport properties of suspensions. With this article, we introduce an¬ 
other method - the renormalized Clausius-Mossotti approximation. Neither of the above 
approximations take the two-body hydrodynamic interactions fully into account. It could 
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be verified by a simple analysis of relevant equations on the two-body level. Consequently, 
the strong hydrodynamic interactions of close particles in suspensions are disregarded in 
both approaches. Therefore, to construct a satisfactory method of calculations of transport 
properties of suspensions, which would take the two-body hydrodynamic interactions into 
consideration, remains an open problem of statistical physics. It is worth noting here, that 
to take the two-body hydrodynamic interactions fully into account in the Beenakker-Mazur 
renormalized fluctuation expansion, a resummation up infinite order is needed. An extension 
of the renormalized Clausius-Mossotti approximation to take the two-body hydrodynamic 
interactions fully into consideration can be done, e.g. by modification of the approximation 
(106), adding the two-body contributions. This type of extension is natural, because it goes 
along the line of a usual, systematic generalization of the Clausius-Mossotti approximation 
[41]. In order to fully grasp the two-body effect in the Beenakker-Mazur method, one needs 
to consider all orders of the fluctuation expansion. It is an important difference between the 
renormalized Clausius-Mossotti approximation and the Beenakker-Mazur method. 


There is another intriguing point from the perspective of statistical physics. The renor¬ 
malized Clausius-Mossotti approximation and the Beenakker-Mazur second order approach 
are similar, because an important element of both methods is an effective propagator. In 
the Beenakker-Mazur method, its role is played by the quantity which depends on 
the volume fraction 0, but does not depend on the distribution of particles, e.g. the two- 
body correlation function [16]. Therefore, this propagator is the same for hard-sphere 
suspension and for suspension of charged particles in equilibrium. On the other hand, the ef¬ 
fective propagator Ge// in the ring expansion, on which the renormalized Clausius-Mossotti 
approximation is constructed, depends both, on the hydrodynamic interactions, and the dis¬ 
tribution of particles. This confrontation gives rise to the question concerning sensitivity of 
both methods to the change in structure of suspension - when, for example, an electrostatic 
interparticle repulsion increases in a suspension and the volume fraction remains unchanged. 
This question in the case of the Beenakker-Mazur method has been answered in some situ¬ 
ations: the Beenakker-Mazur method is rather insensitive to the change in the structure of 
suspension [21, 58]. We are going to address the above questions in our further work. 
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Appendix A: Wave dependent sedimentation 

The wave dependent sedimentation coefficient H (g) describes response of qniescent sns- 
pension, vq = 0 , cf. eqnations ( 11 ), nnder action of a sinnsoidal external force density, 

e.g. 

fext (r) = /oe^Reexp (-fge^ • r) exp (-7712 :|). (Al) 

The above external force density corresponds to the plane wave in the direction of the wave 
vector q = ge^. For a moment, we also introdnce the damping factor exp (— 77 1;2|) and later 
we will take the limit 77 —)■ O’*". The above external force is translationally invariant in x and y 
directions, fext (r) = e^/ext {z). For homogeneous and isotropic suspension that property in- 
dnces the same form for the average velocity held, i.e. (v (r)) = {v {z)). Incompressibility 
condition for the average velocity held (v (r)) = {v {z)) implies {v ( 2 :)) =const. The force 
given by equation (Al) has also the mirror symmetry in 2 : direction, /ext {z) = —fext {—z), 
which, for homogeneous and isotropic snspension, indnces the same symmetry for the veloc¬ 
ity held {v ( 2 :)) = — {v {—z)). Along with the incompressibility condition {v ( 2 :)) =const, we 
obtain, that the average velocity held vanishes in the whole suspension, (v(r)) = 0 . 

This property of the zero net hnx, (v(r)) = 0, simplify the eqnation (26). In this 
sitnation, the npper component of the vector [(U), (f)] in the eqnation (26) is given by 

(U (R, r)) = J d^R'd^r'PuT'^'^ (R, r; R', r') P^iext (r'), (A 2 ) 

were projector Pjj (with its transposition Pf)’), by dehnition projects on the npper half of 
the donble vector [(U), (f)], as follows from the formula (38). From now on, we will take 
the limit 77 —)■ O’*", in which the zero net hux condition, (v (r)) = 0, remains. 
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The external force density given by equation (Al) is torque-free, because the torque acting 
on the particle at position R vanishes, f d^rO (|r — R| — a) (r — R) x fext (r) = 0, for any 
position R. The Heaviside function 6^ (|r — R| — a) used here, vanishes outside the particle 
centered at R and is equal to 1 inside the volume of the particle. It is very simple to show, 
that in the case of the torque free external force density, action of M> operator given by 
equation (9), for such external force density for hard spheres, can be written as 



where Fext (R) is the total force acting on a particle centered at position R, 

Fext (R) = j dhe (|r - R| - a) fext (r) . (A4) 

The equation (A3) is a straightforward consequence of the form of M> given in the reference 
[30]. Similar holds for the single particle Mq operator, 

Ui (Ri; r) = [ d^r'Mo (Ri, r, r') Fext (Ri) • (AS) 

Therefore, action of the operators M> and Mq on the force density fext in equations (10), 
can be replaced using (A3) and (AS). It has the following implication in the equation (A2), 

(U (R, v)) = j^[ d^R' [ (R, r; R', r') P^Fext (R') • (A6) 

^7T (X J J 

In the problem of sedimentation, we consider the average translational velocity of particles, 
(V (R)) = 5 (R — Rj) Vj (X)^ For hard spheres, (V (R)) is related to the average 

velocity held of particles (U (R, r)) , by the formula 

(V (R)> d'r (U (R, r)), (A7) 

|7ra-^ J 

which follows from the relation (2). 

The expression (A6), after passing from the particle velocity held (U(R, r)) for hard 
spheres, to the average translation velocity (V (R)) , yields 

(V (R)) = j d^R'Y (R, R) Fext (R'), (AS) 

where 

Y (R, R') =-^ f d^r 


d^r'PuT" (R, r; R', r') P; 


(A9) 
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For homogeneous suspension, the above kernel Y (R, R') depends only on the difference of 
positions, therefore we can introduce 

Y (R - R') = Y (R, R'). (AlO) 

Moreover, isotropy implies, that Y (R) is of the form uq (R) l+y 2 {R) RR. This form of 
Y (R), along with simple algebraic manipulations of equations (AS), leads to the following 
conclusion. For the external force given by a plane wave, 

Fext (R) = FoqReexp (-iq ■ R), (All) 

the average translational velocity has also the plane-wave form, 

(V (R)) = R (q) q Re exp (-fqR), (A12) 

where the coefficient V (g) is given by the formula 

V (q) = H (q) yoFo, 

with hydrodynamic function given by equation (36). The stokes coefficient po = l/{6nria). 

Appendix B: Calculations within renormalized Clausius-Mossotti approximation 

In this appendix we give some details of our calculations within the renormalized Clausius- 
Mossotti approximation. It demands to solve the set of equations (106), (105), (27), (78), 
(75) and (77) for the quantities B, Gg//, T*’’’’, T, Tog. Each of those quantities is 

a fourfold 3x3 matrix, as it may be inferred for example from the formula (22) for the T 
operator. Therefore the T operator written with all variables and indexes is represented by 
T(R, r, R', r%Q,ti/^Q,/. R has two indexes u,u' = U,P, which denotes upper or lower part of 
the double vector, and another two Cartesian indexes a, a' = 1, 2, 3. In our calculations we 
represent those quantities as multipole hydrodynamic matrices [59]. In the reference [30] the 
reader can End how to introduce a multipole picture for the forces and velocities and how 
to represent the hydrodynamic matrices M and G in the multipole picture. The notation 
used in the reference [30] is adopted also in this article, therefore we do not repeat that 
material. In the same way as in the reference [30] we introduce the multipole picture of 
the hydrodynamic matrices B, Ge//, T*’’’’, T and Toff. In the multipole picture, all 
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of the above quantities become infinite dimensional hydrodynamic multipole matrices, e.g. 
T (R, r; R', [T Therefore the variables r, r' and Cartesian 

indexes are transformed into multipole numbers I, m and a having the following range: I = 
1,2,... ,oo, m = —I, + whereas a = 0,1, 2. An important role in onr calculations 

is played by the homogeneity of the system, because it implies that the matrices depend on 
the difference of positions only, for example [T (R - = [T (R, 

We also use isotropy of the system, which allow to calculate the multipole matrix T (R) for 
any vector R when the T for direction R = is known. 

Using also the Fonrier space, with the Fonrier transformation given by 

f (k) = J d^R exp[-ikR]T(K), (Bl) 

with similar dehnition for other multipole hydrodynamic matrices, we represent the set of 
eqnations (106), (105), (27), (78), (75) and (77) as follows 


where 


T(k) = R + Toff(k), 

(B2) 

f (k) = r’'’'(k) (i - G(k)r"(k))”^, 

(B3) 

Gefr(k) = G(k) + G(k)f(k)G(k), 

(B4) 

f‘^^M(k) ^ B, 

(B5) 

’^’■’'(k) = T^6M(k) (l - [7S](k)f^^M(k))'', 

(B6) 

B = riiM + J d^R Toff(-R)G'(R)M, 

(B7) 

{hG^]{k) = j d^R ea;pHkR]h(R)Geff(R). 

(B8) 


In the above equations there appear snperpositions of the mnltipole hydrodynamic matrices, 
inversion of matrices, and the matrices appear both in the positional and the Fonrier space. 

We created numerical code to solve these equations of the renormalized Clausius-Mossotti 
approximation. Three aspects appear here. 

First, in our numerical calculations, we truncate the hydrodynamic matrices, e.g. T (k) 

L J ulmcr,u'l'm'a' 

The trnncation is characterized by L, which is the highest multipole used in the calculations 
- we consider matrix elements with 1,1' < L only. The calculations were performed for 
different parameters trnncation, L = 4,..., 10. The dependence of the effective viscosity 
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coefficient T]es/v on the function 1/ (Llog^L) of the truncation parameter L is presented in 
the hg. 6. The hgure shows, that even for the highest truncation parameter L = 10 (which 
corresponds to 1/ (Llog^L) 0.008) the effective viscosity coefficient is still sensitive to 
the truncation L. Therefore, we extrapolate the coefficient up to L ^ cxd. The extrapolated 
value of the coefficient is given by intersection of a straight line, passing through the points 
corresponding to L = 9 and L = 10 in the £g. 6, with the vertical axis. Similar procedure 
of the extrapolation is carried in the case of the other transport characteristics. 

The second numerical aspect is related to a discretization of distance R for hydrodynamic 
matrices, e.g. T Points R = ..., with .^o = 0 and = ■Ci exp [a (i — 1)] were 

considered. The parameter a was determined from the assumption, that the first and the 
last section are equal, i.e., —Therefore, .^i and N determine a set of points, in 

which the hydrodynamic matrices were considered in the code. Other values - if needed - were 
calculated by interpolation or from the asymptotic expansion (e.g. for R larger than It 
is worth noting that the exponential mesh is convenient to calculate the three dimensional 
Fourier transform of the hydrodynamic matrices. In our calculations the three dimensional 
Fourier transform of hydrodynamic matrices in the formula (Bl) was first reduced to the 
one dimensional Hankel transform [60] in a similarity to a dielectric system [61]. Then the 
exponential mesh is used to perform Hankel transform with the use of numerical procedures 
for fast Fourier transform. We performed calculations for ^i/ (2a) = 1/2,1/3,1/4,1/5 and 
N = 512,1024, 2048,4096 respectively. Larger N correspond to a denser mesh. We found 
the mesh characterized by iV = 4096 sufficient and it is used to obtain the results presented 
in this article. 

The third aspect of the numerical calculations is related to the fact, that the system of 
equations was solved iteratively. We observed, that up to the volume fractions 0 ~ 45%, 
after a few iterations, difference between hydrodynamic functions in subsequent iterations 
decays as it happens in a geometric series. It may be written as follows 

sup iHi (q) - Hi_i {q) \ ^ sup \Hi+i (q) - Hi {q)\ A, (B9) 

g g 

where Hi (q) denotes the hydrodynamic function after the i-th iteration and symbol sup^ / (q) 
stands for maximal value of a function /. The highest observed value of A 0.7. The 
iteration procedure was stopped when the following condition sup^ \Hi (q) — Hi_i (g)! < 10“"^ 
was satisfied. 
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It is worth mentioning, that the computer time and memory to solve the equations 
iteratively is comparable with the calculations, which we performed in the case of the revised 
Beenakker-Mazur method [30]. The numerical results within the renormalized Clausius- 
Mossotti approximation presented in this article were calculated with a use of a desktop 
computer within one day. 
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FIG. 6. The relative effective viscosity coefficient 'Heff/'n as a function of multipole truncation L 
obtained by the renormalized Clausius-Mossotti approximation for suspension of volume fraction 
(/. = 0.45. 
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